Phase-Field Modeling of Nonlinear Material 
Behavior 

Y.-P. Pellegrini, C. Denoual and L. Truskinovsky 



Abstract Materials that undergo internal transformations are usually described in 
solid mechanics by multi-well energy functions that account for both elastic and 
transformational behavior. In order to separate the two effects, physicists use in- 
stead phase-field-type theories where conventional linear elastic strain is quadrati- 
cally coupled to an additional field that describes the evolution of the reference state 
and solely accounts for nonlinearity. In this paper we propose a systematic method 
allowing one to split the nonconvex energy into harmonic and nonharmonic parts 
and to convert a nonconvex mechanical problem into a partially linearized phase- 
field problem. The main ideas are illustrated using the simplest framework of the 
Peierls-Nabarro dislocation model. 



1 Introduction 

Nonconvex energy potentials are used in solid mechanics for the modeling of 
martensitic transformations [9], plasticity [1] and fracture [25]. Parts of the resulting 
energy landscapes correspond to sufficiently smooth deformations preserving the lo- 
cally affine structure of the lattice environment of each atom. Other parts represent 
highly distorted atomic arrangements associated with either loss or reacquisition of 
nearest neighbors. While deformations of the first type can (often) be described by 
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the conventional strain tensor of (linear) elasticity theory, a representation of the 
deformations of the second type requires introducing additional internal variables 
accounting for deviations from the local affinity of the stressed atomic configura- 
tions. In particular, these supplementary variables describe the evolution of the local 
reference state (LRS) from which the elastic deformations are measured [3, 1 1, 26]. 
The main difference between the elastic strains and these supplementary internal 
variables is that the dynamics of the former is typically inertial, while that of the 
latter is usually overdamped. Sometimes the nonelastic variables can be minimized 
out as in the case of deformational plasticity (e.g., [2]). In this paper we deal instead 
with situations where the internal variables have to revealed rather than hidden. 

We assume that the coarse-grained nonconvex energy density /(e) is known ei- 
ther from extrapolations of experimental measurements or from ab-initio calcula- 
tions involving atomic homogeneity constraints. We suppose that the argument e 
of this function, that represents a coarse-grained strain, is small and can be addi- 
tively split into the linear elastic part e, and a phase-field part rj that accounts for the 
nonelastic evolution of the LRS. Our next assumption is that / can be represented 
as a sum of two terms: the elastic energy f e , which depends on e = e — 77 and the 
phase-field energy g, which depends on rj. We interpret /(e) as the outcome of adi- 
abatic elimination of the variable rj and consider the inverse problem of recovering 
the phase-field energy #(17) from the function /(e) under the assumption that the 
function /.(e) is quadratic. The problem of the identification of g(Tj) reduces to a 
problem of optimization and the relation between the 'optimally' related functions 
/(e) andg(rj ) is studied in some prototypical cases. If, in contrast, the function g(rj) 
is chosen independently, the corresponding function /(e) is typically non-smooth 
and non single-valued, e.g. [6]. 

To motivate the need for the phase-field variables we consider in full detail a spe- 
cific physical example. It deals with the mixed, discrete-continuum representation of 
a dislocation core [12, 16]. More specifically, we develop a modified version of the 
classical Peierls-Nabarro (PN) model that accounts for a finite thickness of the slip 
region. In this problem the coarse-grained description of the slip zone is provided 
by the so-called /-potential [5, 27]. The phase field represents an "atomically sharp" 
slip and the part of the interaction potential related to g gives rise to the slip-related 
pull-back force [7, 16, 23]. Our general method of recovering the expression for this 
force represents an extension of Rice's transform, which was first introduced in the 
context of a dislocation nucleation problem [19]. 

In this paper only the simplest scalar problem in a one-dimensional setting is con- 
sidered. The slightly more general question of extracting from the coarse-grained 
energy a convex (instead of quadratic) component will be examined elsewhere [17]. 



2 Surface Problem 

We begin with the special case when the phase field is localized on a surface. In 
problems involving fracture or slip it often proves convenient to represent the energy 
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of a body as the sum of a bulk term depending on strain gradients and a surface 
term penalizing displacement discontinuities. The bulk term is usually modeled by 
linear elasticity. The modeling of the surface energy is less straightforward [6, 25]. 
For instance, the models will be different depending on whether the location of the 
discontinuities is known a priori or not. 

In a ID setting with a known fracture set the equilibrium problem reduces to 
minimizing the following energy functional 



Here f e {e) = (E /2)e 2 , E > is the elastic modulus and a is a coarse-graining 
length scale that typically exceeds several atomic sizes. The set T a in (1) represents 
discontinuity points resolved at scale a and 8(x) = [u] a (x] is the corresponding dis- 
placement discontinuity. The surface energy f a (8) is then an effective interaction 
over the distance a; in particular, the shear-related component of f a (8) coincides 
the /-potential mentioned in the Introduction. 

In the case when the fracture set is unknown the surface energy has to be chosen 
differently. The reason is that in this model the displacement discontinuity at scale 
a does not represent the microscopic slip between neighboring atomic planes, and 
therefore the difference between elastic deformation and inelastic slip has to be yet 
resolved at this scale [19]. More precisely, linear elasticity, which has nothing to do 
with slip and which is already accounted for in the bulk term, has not been excluded 
from f a (8). The identification of the surface energy in (1) with f a (8), which is 
quadratic at the origin, leads in a free discontinuity problem to a degenerate solution 
with infinitely many infinitely small discontinuities [6]. 

To remove linear elasticity from the surface term, one should replace the coarse- 
grained discontinuity [u] a by the atomically sharp slip rj(x) = [u](x) that does not 
depend on a. The energy (1) is then rewritten as 



where now F is the set of discontinuity points corresponding to a = 0. The problem 
is to find the relation between the function f a (8), representing an empirical input, 
and the unknown function g(rj). 

To define g(rj) we divide the total slip 8 into an elastic part, ae, where e is an 
equivalent elastic strain, and an inelastic part 77. The function g(r\) is defined by the 
condition that f a (8) is arelaxation of the energy af e (e) + g(rj) under the condition 
that ae + 77 = 8, namely: 




a 



(1) 




(2) 




(3) 
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If the energy f a {8) is a single-well function and the infimum is unique, the func- 
tion g is completely defined. If f a (8) is periodic as in the case of dislocations, 
in order to have a uniquely-defined g(ri), we need to replace in definition (3) the 
global minimization by a properly-defined local minimization denoted hereafter by 
'infioc' (minimization over 77 starting from the minimum of f a closest to 8.) 

In what follows, our task will be to reverse definition (3) and to recover the 
nonequilibrium energy g(rj) from f a (8). What allows us to proceed is the specific 
(harmonic) structure of the elastic part of the energy. 

We observe that the function g must satisfy the following necessary condition 

(E/a)(8-r 1 )=g'(r 1 ). (4) 

Moreover, differentiation of (3) wrt. 8 gives 

f a {8) = {E/a){8-n). (5) 

These two equations allow one to represent g'(rj) in the following parametric form 
[7, 8, 19] 

(r],g'(T 1 )) = {8-(a/E)f>(8),f>(S)). (6) 
The parametric representation for g(Tj) then reads 

= (S {a/E)f a {8)J a {8) [a/(2E)} [f^S)] 2 ) . (7) 

Since for nonconvex f a (8) this representation may lead to a multivalued function 
g(rj) formula (7) must be supplemented by an additional branch selection proce- 
dure. To illustrate the mapping f(8) —> g(rj) given by (7) and the selection of a 




Fig. 1 Parametric transforms (6), (7) applied to a Lennard- Jones potential f(S) = 8~ 12 — 8~ 6 . a) 
f(S) (dashed), and the two branches of g(rj) (solid), where 'p' ('u') labels the physical (unphysi- 
cal) branch; b) f'(S) and the 'p' branch of g'(T]); c) T](5). 



physical branch we consider a Lennard- Jones potential /„, with a = 1 and assume 
that E = f"(So), where 80 is the only minimum of / (see Fig. 1). Notice that the 
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resulting function g'(ri) has an infinite slope at 77 = 80 and that for 77 > 8q we must 
haveg(T?) oc (77 -<5o) 3//2 . 

The removal of the linear elastic part of the energy becomes important in PN- 
type modeling of dislocations. Consider, for instance, a straight screw dislocation in 
an isotropic linear-elastic body and assume that the sharp discontinuity plane, y = 0, 
lies between the two effective gliding surfaces located at y = ±a/2. To account 
for the finite thickness of the core region a we need to modify the classical PN 
model [12]. According to our interpretation the linear elastic stress outside the slip 
region (—a/2, a/2) must be balanced by the coarse-grained pull-back stress that is 
resolved at the spatial scale a. We therefore interpret the pull-back stress at this scale 
as f^(8(x)), where f a is the y-potential, a periodic function with period b and with 
fa (0) = 0- The expression for the linear stress outside the slip region is derived in the 
Appendix. With these considerations in mind we obtain for the unknown function 
tj (x) representing a mathematical slip at y = the following system of equations 

~ xr^' arctan 2o^~7) + ° a(?c) = fa - (8) 

8{x) = (a/ i i)f' a (8(x))+T]{x), (9) 

where (7„ is the resolved applied stress at scale a. If we match the linear elastic 
behavior at rj = with that in the bulk regions we obtain that /i = af"(0). Using in 
this relation the physical shear modulus and the value of f" (0) from the y-potential 
provides a rough estimate for a, the effective interaction range. 

We notice that parameter a enters both equations (8,9), which makes this system 
different from the one studied in [16, 19]. The ideas behind our nonlocal extension 
of the PN model are also different from that of Ref. [13] where a nonlocal kernel 
was introduced empirically as part of the pull-back stress, and the usual l/(x — xf) 
kernel was used for the bulk stress. 

To bring the system (8,9) into the framework of phase-field models, we identify 
the effective pull-back force f' a (5(tj)) with g'{r\) and rewrite Eq. (8) as 

^J^c'K a (x-x')7]'(x') + a a (x) =*'(tj). (10) 

where K a (x) = — (2/7z;a)arctan(a/2jc). It is now easy to see that g'(rj) enjoys the 
parametric representation 

(n,tf(Ti)) = (8-~fa{8),f a (8)), (ID 

where we recognize the mapping (6) (see also [16, 19, 23]). To make the link with 
the classical PN model one needs to consider the limit a — > 0. By computing T7 in 
terms of 8 and expanding (10) in powers of a, we obtain to order 0(a) the following 
'gradient' extension of the PN model 
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| + Jck'^ + A5"(x) + a f ,(x) = /(«(*)), (12) 

where A = aji/A. For different weakly or strongly nonlocal generalizations of the 
PN model see [13, 21]. Equation (12) features an effective applied stress that dif- 
fers from (7„(x,0) defined in the Appendix by an 0(a) correction, namely a a {x) = 
a a (x,0) + (a/2) [(l/2)d y a (x,O)-Ko*d x o a (x,O)]. The classical PN model is re- 
trieved by letting a = 0. 



3 Bulk Problem 



Now let us place the problem in a more general framework. The task is to approxi- 
mate locally the empirical potential /(e) by a quadratic function with an optimally 
chosen reference state rj, and to associate with this state a reference energy g(rj). 
Behind such construction is the assumption that all the nonlinearity of the problem is 
related to the evolution of the reference state. The simplest setting to pose formally 
the problem is the one-dimensional geometrically linearized theory of nonlinear 
elastic bars. 

According to our interpretation the empirical energy is represented as 



/(£) = inf 

7],l0C 



|(e-77) 2 + g(77) 



(13) 



and the problem is to find the intrinsic phase-field function g(f]). Following the 
previous section we write the parametric representation for g'(r\) in the form 

(n,g'(n)) = (e- i ^ 1 ,f(e))- (14) 

The function g(rj) is then given by the mapping 

fa.*fo)) = («-^./(«)-^)- ( 15 > 

The consistency of this procedure requires the parameter E and the function /(e) 
to be related. If we expand the parametric definition of g(rj) near a reference state 
e where /'(eo) = 0, we obtain g(e ) = /(e ), g'(e ) = and g"(e ) = /"(e )/[l - 
/"(eo )/E]. The natural choice E = /"(fio) makes g"(eo) infinite. The behavior of 
the higher derivatives of g(rj) near 77 = eo is found by assuming (without loss of 
generality) that derivatives /^(£b) vanish for k = 3,...,n— 1. The order of the 
asymptotics depends on n > 2, which is the first integer such that /'"'(eo) 7^ 0: 
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Hence g behaves near its minimum as: \g(rj) — g(£o)\ ~ |lj — £o\"^" ~ ■ The generic 
case is « = 3; the case n = 4 corresponds to periodic potential relevant for disloca- 
tions; for / locally harmonic, n = °°. 




-i-i—i 

e T| £(T() 



Fig. 2 Geometrical illustration of the construction defined by Eqs. (14-17). 



Observe now that the function g computed from (14,15), can also be viewed as a 
solution of the following optimization problem: 



g(r]) = sup 



£,loc 



(17) 



which is a natural inverse of (13) (see also [18, 20]). Since the equation 77 = 
e — f'(e)/E may have several solutions e(rj), the representation (17) removes the 
ambiguity by always selecting the upper branch. 




8 



Y.-P. Pellegrini, C. Denoual and L. Truskinovsky 



The working of Eqs. (14-17) with E = f"(eo) is illustrated in Fig. 2. In the do- 
main at the left of £o, where / grows faster than harmonic the desired tangency point 
does not exist. In this case the difference /(e) — j (e — rj ) 2 is maximized at e = — °o. 
This situation takes place in the Lennard- Jones example of Sec. 2 where we have to 
use g(8) = +°° for 8 < 8q (hatched area of Fig. la). 

To handle general multi-well energies, we first introduce the Stillinger- Weber 
mapping £o(e) that links to any state e the local minimum £o of /(e) that would 
be attained from this state by steepest-descent [22]. Next, we modify equations (13) 
and (17) as: 



inf 

7],l0C 



g(~n) = sup 

e.loc 



r( eo ( e ))( e -77) 2 +s(77) 



m-^f"(eo(n))(e-nf 



(18) 
(19) 



Whereas (19) defines g, equation (18) states that knowing / is equivalent to knowing 
g plus the linear-elastic behavior of / near its local minima. 

The precise meaning of the "loc" in Eqs. (18), (19) is as follows. Operationally, 
the minimization in the definition of / is carried out over r\, starting from £o(e), 
the local minimum nearest to e determined by the SW mapping; the corresponding 
elastic modulus is also determined by the starting point. The maximization in the 
definition of g{v\) proceeds along similar lines except that now the relevant elastic 
modulus is determined by the local minimum closest to T], and is fixed during the 
maximization. The optimization is carried out starting from e = rj . 

Figs. 3 illustrate the case of a double-well potential with unequal curvatures of 
the wells. Notice that in contrast to what we saw in Fig. 1 the function 77(e) is now 
bounded. 




Another interesting case is the periodic potential that is used in the description 
of reconstructive phase transitions (e.g., [4]). Consider, for instance, the piecewise- 
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harmonic periodic case shown in Fig. 4 that is often used in analytical studies [13, 
21]. The parametric representation (14) of g is here useless and the definition (19) 
must be used instead. In this extreme case, all elasticity has been removed from g 
and the resulting g(r\) is cone-shaped at its minima (Fig. 4a) as predicted by Eq. 
(16) for n — > oo. The force g(rj) is discontinuous (Fig. 4b) and its extreme values 
provide thresholds for the evolution of 77, whose stepwise character is an artifact 
due to the absence of smooth spinodal regions in /. 

It is also instructive to consider for comparison the case of an unbounded har- 
monic potential /(e) = (Ef/2)(e — So) 2 . From (17) with E = Ef, one deduces that 
gitf) = +°° if ^7 7^ £o an d gill = £o) = 0. This trivial example indicates that in a 
purely linear-elastic model, the reference state does not have to evolve. 



4 Concluding Remarks 

The goal of this paper was to reveal in the simplest setting the variational nature of 
the generalized Rice transform. The problem consists in splitting a coarse-grained 
lattice potential /, describing the overall deformation of a sufficiently large number 
of atoms, into a (quasi) convex elastic potential and an inelastic potential g dealing 
with structural rearrangements. Here the potential / is assumed to be measurable 
by molecular statics along a prescribed deformation path relevant to the material 
transformation in question. For simplicity, the elastic potential is assumed in this 
paper to be a standard quadratic function of the macroscopic strain. The inelastic 
potential g must be a function of the phase-field variable 77, whose identification 
represents an important part of the problem. While our precise construction solving 
the above problem is presented in the static setting (see Eq. (19)), the motivation for 
the splitting concerns, first of all, dynamical applications (e.g., [7]). Thus we assume 
that material displacement u associated to the strain e = d x ii evolves inertially almost 
without damping (standard elastodynamics), while the dynamics of the phase-field 
variable is overdamped. More precisely we assume that the relaxation of the variable 
77 follows the time-dependent Ginzburg-Landau (TDGL) equation. By means of an 
empirical 'viscosity' parameter v we can write the evolution equation in the form 



where we have omitted for simplicity the conventional gradient-penalizing terms 
(e.g., [7, 24]). In the static setting the above equation reduces to our basic Equ. (18). 
The definitive advantage of separating the wave motion from an overdamped TDGL 
relaxation is the possibility to attribute effective damping only to large atomic dis- 
placements. Our preliminary investigations [17] indicate that extending the varia- 
tional set-up presented in this paper to higher dimensions and generalizing it in the 
direction of extracting (quasi) convex, rather than merely quadratic elastic compo- 
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nents, is feasible. These issues will be addressed systematically in a separate publi- 
cation. 



Appendix 

The following computations are largely based on the Eshelby's arguments pre- 
sented in [10]. Consider a Volterra screw dislocation with zero-width core and with 
Burgers vector b. The displacement u z (x,y) has the form 

b b y b 

Uz(x,y) = — Arg(x + iy) = — arctan - + - sign(y) d(-x), (20) 
2n 2n x 2 

where 9 is the Heaviside function, and where the indeterminacy in the discontinuity 
of u z is resolved by specifying the glide plane (y = 0). The distributional part in 
the r.h.s. of Eq. (20), usually omitted in the literature (e.g. [12]), is crucial to the 
present derivation because it represents the irreversible atomic displacements on the 
plane y = 0. We introduce the eigendistortion, as the part of the dislocation- 
induced distortion J3,y = uij, that is not linear-elastic. For our dislocation, its only 
non-zero component is j3*-(x,y) = bQ{ — x)8o{y), where 8d is the Dirac distribution 
[14]. The linear-elastic distortion, J3?-, is defined through the additive decomposition 
of the total distortion j3,y, namely j3?- = J3y — /3 ; *- [14]. The elastic strains are e ;/ - = 
sym/3/y. By using the identity [arctan (1/x)]' = k8d{x) — 1/(1 +x 2 ), we obtain that 
the distributional parts in J3,y and j3 ; *- mutually cancel out giving the standard result 
[12] ^ ^ 

?xz{x,y) = ~r~ 2 , 2 ' e yz{x,y) ~ ~T~ 2 T 2 ' ( 21 ^ 
47Tx z +y z 4^x +y 



The stress induced by the eigenstrain is then a*(x,y) = 2/ie,v(x,y), where i = x,y. 
In the presence of an applied shear stress [15] O a = Oayz, Eq. (20) becomes 

u : (x,y) = - [ } dy'<j a (x 7 y') + ^- arctan^ + ^-sign{y)e(-x), (22) 
jl Jo 2n x 2 

The total stress is then a = a* + a a and e = g/2/j.. Now, the key step consists in 
averaging the stress over the layer of width a containing the glide plane. Introduce: 
cjjj (x) = i f^"/2 dy °ij ( x ^ y) ■ From (22), the x component of the total relative atomic 
lattice displacement between the atomic planes at y = ±a/2 reads: 

■ffo0(-x). (23) 



a 



8{x) = u z (x, +a/2)— u z (x, —a/2) 



_ , , lib ( a \ 

<y a {x)-\ arctan — 

Ka \2xJ 



Furthermore on account of (21) the average shear stress O yz (x) in the layer is: 

lib 1 c a l^- x lib i cl \ 

o yz {x) = <r fl W + --y d y?T? = a„(x) + -arctan (-) . (24) 
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Comparison of (24) and(23) shows that: 



S(x) = (a/n)a yz (x)+bd(-x). 



(25) 



Consider next an Eshelby screw dislocation with an extended core described by 
a continuous function t](x). The distortion j3* becomes: fi*(x,y) = 5d (y) 77 (jt) = 
— Soiy) J x 6x7]' (x); the Volterra dislocation corresponds to the limiting case 17 (x) = 
b6(—x). Displacements, strains and stresses are obtained by convolution using 
dj3* (x,y) = — 8D(y)ri'(x)dx [10] as elementary distortions. The analogs of Eqs. 
(24,25) are: 



Equation (27), which we use in the paper, shows the relation between the coarse- 
grained displacement 8, and the discontinuity 77. 
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